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1.  Introduction 


The  classical  trajectory  problem  in  a  uniform  gravitational  field  is  well  known,  if  there  are  no  drag 
forces  on  the  projectile.  The  trajectory  traces  out  a  parabolic  path  in  the  atmosphere  over  time, 
while  the  horizontal  component  of  velocity  remains  undiminished.  In  the  real  world,  however, 
the  projectile  is  subject  to  drag  forces,  which,  at  high  Reynolds  number,  are  estimated  to  be 
proportional  to  the  square  of  the  velocity  ( i.e .,  aerodynamic  drag).  This  drag  force  complicates 
the  solution  of  the  trajectory. 

Prior  studies  attempting  to  analytically  address  the  issue  have  been  limited  to  a  variety  of 
approximate  methods,  e.g.,  (7).  Alternately,  others  (2)  have  modified  the  assumption  of 
aerodynamic  drag  in  favor  of  a  (less  physically  based)  constant-deceleration  approximation. 

Along  similar  lines,  some (3)  redefine  the  drag  as  viscous  (i.e.,  proportional  to  velocity  V ),  rather 
than  aerodynamic  (proportional  to  V2).  Still  others  (4)  develop  separate  equations  for 
pure-vertical  and  pure-horizontal  flight,  under  the  influence  of  aerodynamic  drag,  and  then 
assume  (by  way  of  approximation)  that  they  can  be  respectively  applied  to  the  vertical  and 
horizontal  components  of  an  arbitrary  trajectory.  With  the  advent  of  a  pervasive  computing 
environment,  numerical  approaches  can  be  adopted  to  calculate  trajectory  solutions.  However, 
numerical  solutions,  which  apply  only  to  a  particular  set  of  initial  conditions,  fail  to  provide  the 
functional  understanding  of  trajectory  behavior  afforded  by  an  analytical  solution. 

A  prior  analytical  study  was  conducted  (5)  to  determine  the  distance  a  hypervelocity  projectile 
would  travel  in  air.  In  that  study,  special  cases  of  the  more  general  problem  were  analytically 
solved,  such  as  vertical  launch  and/or  small- angle  (i.e.,  near  horizontal)  trajectory.  In  the  current 
effort,  the  more  general  case  of  arbitrary  launch  angle  and  large  changes  in  the  trajectory  are 
considered.  The  projectile  is  ballistic,  by  which  we  mean  capable  of  neither  propulsion  nor  lift. 
The  model  is  based  on  Newton’s  2nd  Law  and  the  sole  forces  acting  on  the  projectile  are  drag  and 
gravity.  Erosion,  ablation,  and  strength  of  the  projectile  are  not  considered.  This  study  is  useful 
for  experimental-range  safety  considerations,  as  well  as  fragment  lethality  and  fragment  recovery. 

The  problems  of  range  safety  and  danger  areas  have  been  examined  in  the  past  (6),  with 
aerodynamic  drag  being  applied  to  various  ballistic  projectiles.  However,  such  prior  studies  have 
not  put  forth  analytical  relations,  such  as  those  pursued  in  the  current  work. 
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2.  Flight  Retardation  Equations  (Cartesian  Coordinates) 


Consider  a  projectile  of  mass  m  launched  at  an  initial  velocity  V  =  V0,  at  an  angle  o0  with  respect 
to  the  horizon,  subject  to  the  initial  (time  t  —  0)  conditions  that  the  horizontal  position  x  —  0,  and 
the  vertical  position  y  takes  on  a  value  of  H,  representing  the  height  above  ground  of  the  launch. 

The  velocity  components  at  any  given  moment  are  given  as 

vx  —  x  —  V  cos  a  (1) 


and 


Vy  —  y  —  V  sin  a 


(2) 


where  a  is  the  time-dependent  trajectory  angle  with  respect  to  the  horizon,  and  the  overdot 
denotes  time  differentiation.  The  aerodynamic  drag  force  on  the  projectile  is  pACdV2 / 2,  where 
p  is  the  density  of  air  at  sea  level  (1.293  kg/m3),  g  is  the  acceleration  due  to  gravity  (9.806  m/s2), 
Cd  is  the  drag  coefficient,  and  A  is  the  effective  cross-sectional  area  (m2).  The  governing 
equations,  accounting  for  the  effects  of  aerodynamic  drag  and  gravity  and  neglecting  ablation, 
according  to  Newton’s  2nd  Law,  are: 


pACdV 2 

rnx  = - - -  cos  a 


pACd. 
— - — x 


(3) 


my  = 


pACdV 2  .  pACd.  r 


sm  a  —  mg  = 


ijyx2  +  y2  —  mg  . 


By  lumping  the  term,  B  =  pACd/2m ,  we  may  restate  the  governing  set  of  equations  (in  Cartesian 
coordinates)  as 

x  =  — BxJx 2  +  y 2  (5) 


V  =  -BijJx2  +  y2  -  g  . 


The  grouping  within  B,  given  by  m/A,  represents  the  areal  density  along  the  flight  axis  of  the 
projectile,  and  may  be  alternately  expressed  as  ppLp,  where  pp  is  the  projectile  density  and  Lp  is 
the  characteristic  length  of  the  projectile.  Thus,  an  alternative  expression  for  B  is  given  by 
B  =  pCd/2ppLp.  Note,  that  in  reference  4  equation  5  is  solved  assuming  y  —  0,  and  equation  6 
is  solved  assuming  x  —  0,  resulting  in  greatly  simplified,  but  approximate,  solutions. 
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3.  Flight  Retardation  Equations  (Curvilinear  Coordinates) 


Define  direction  s  as  the  coordinate  along  the  flight  trajectory,  and  n  as  the  direction  normal  to  the 
current  trajectory.  Because  the  curvilinear  coordinate  system  follows  the  trajectory,  the 
coordinate  n  will  always  be  identically  zero,  even  as  there  are  forces  bending  the  coordinate 
system  in  the  n  direction.  For  example,  there  is  a  component  of  the  gravitational  force  in  both  the 
s  and  n  directions.  Additionally,  there  is  aerodynamic  drag  along  the  axial  s  direction.  The 
governing  equations,  in  curvilinear  coordinates,  are  therefore 

s  —  —gsma  —  Bs2  (7) 

n  —  g  cos  a  .  (8) 


The  tendency  of  the  trajectory  (and  the  coordinate  system)  to  bend  by  way  of  equation  8  may  be 
related  to  changes  in  the  trajectory  angle  a.  From  a  velocity  diagram,  figure  1, 


—  tan  da 


n  dt 

- — da 

s 


(9) 
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Employ  a  =  da/dt  and  rearrange: 


(10) 


— ds  =  n 


Significantly,  this  equation  relates  lateral  forces  on  the  projectile  to  the  speed  and  changes  in  the 
trajectory  of  the  projectile.  Substitute  equation  8  to  eliminate  h: 

as  =  —  g  cos  a  .  (11) 

Equation  1 1  essentially  replaces  the  n  relation  of  equation  8  as  the  second  governing  equation,  in 
addition  to  equation  7.  We  thus  have  two  curvilinear  governing  equations,  in  terms  of  s  and  a 
and  their  derivatives. 


4.  Relating  s  and  a  Without  Curvilinear  Coordinates 


Equation  1 1  can  be  alternately  derived  without  the  use  of  curvilinear  coordinates,  as  follows. 
Multiply  equation  6  by  x: 

x(y  +  g)  =  —Bxy\J  x2  +  y2  (12) 

and  equation  5  by  y: 

ijx  =  — Bxi)\jx 2  +  y2  .  (13) 

Subtract  equation  13  from  equation  12: 

xy  —  yx  +  gx  =  0  .  (14) 


Rearrange: 


Set  this  result  aside  for  a  moment. 


The  trajectory  angle  a  is  defined  as 


Differentiate  with  respect  to  time: 


xy  —  yx  g 

x2  x 


tana  = 


y 

x 


9  xy  —  yx 

a  sec  a  = - - — 

x 2 


(15) 


(16) 


(17) 
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One  can  use  the  result  of  equation  15  to  eliminate  the  right-hand  side  of  equation  17: 


2  9  9  sec  a 

a  sec  a  =  —  = - 

x  s 


where  the  geometric  pathlength  relation  ds/dx  =  sec  a  has  been  used.  This  result  becomes 


g  cos  a 


a 


which  is  identical  to  the  result  of  equation  11. 


(18) 


(19) 


5.  Second-order  Governing  Relation  in  One  Variable 


Take  equation  19  and  differentiate  with  respect  to  time: 


g  cos  a  a 

s  =  g  sin  a  H - „ - 

cr 


Equate  this  with  equation  7  to  eliminate  s: 


(20) 


2  g  cos  a  a 

-g  sm  a  —  Bs  —  r/sm  «  -* - -= - 

a1 


(21) 


Employ  equation  19  to  eliminate  s  by  substitution,  which  leaves  the  equation  wholly  in  terms  of 
trajectory  angle  a  and  its  derivatives: 


a  +  2  tan  a  a2  +  Bg  cos  a  =  0 


(22) 


This  second-order  governing  relation  can  also  be  rewritten  as 

d2 

—■(tan  a)  +  Bg  sec  a  =  0  .  (23) 

dt 2 

We  see  that,  when  there  is  no  drag  ( i.e .,  when  B  =  0),  equation  23  yields  a  solution  in  which 
tan  a  is  linear  in  time.  This  corresponds  to  the  classical  parabolic  solution  in  which  x  is  constant 
and  y  is  linear  in  time  (since  tan  a  =  y/x). 
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6.  Integrating  the  Second-order  Governing  Equation 


For  convenience,  assign  the  substitution  w  =  tan  a.  Then,  the  first  term  of  equation  23  is  merely 
it),  which  can  alternately  be  expressed  as  w  dw/div.  In  that  case,  equation  23  may  be  reexpressed 
as 

w  dw  =  —Bg  sec  a  dw  .  (24) 

However,  knowing  the  definition  of  w,  it  follows  that 

dw  =  sec2  a  da  ,  (25) 


(26) 


(27) 


where  C  is  an  integration  constant.  We  choose  to  treat  the  large  bracketed  term  on  the  right-hand 
side  as  an  intermediate  variable,  call  it  u: 


u(a )  =  C  —  tan  a  sec  a  —  ln(tan  a  +  sec  a)  ,  (28) 

such  that  w  =  d/dt (tan  a)  =  —y/Bgy/u.  The  minus  sign  of  the  square  root  has  been  taken,  as 
we  know  that  tan  a  must  be  decreasing  throughout  the  trajectory.  We  further  employ  the  chain 
rule  to  express  vj  as  a  dw/ da  =  a  sec2  a.  Thus,  we  can  present  equation  27  as 


a  =  -^Bg.^-  .  (29) 

v  secz  a 

Since  u  is  a  function  of  a  alone,  we  note  that  equation  29  presents  the  time  rate  of  change  of  the 
trajectory  angle  (d)  in  terms  of  the  trajectory  angle  itself  (a). 


We  may  solve  for  the  integration  constant  C  from  the  launch  conditions.  Equation  1 1  tells  us 


that,  at  launch, 


g  cos  a0 


(30) 
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This  may  be  employed  in  equation  29,  at  the  t  —  0  condition,  to  obtain  u(a0)  and  thus  C  as 


g  sec  a0  . 

C  =  — - b  tan  a0  sec  ao  +  m(tan  a0  +  sec  ao)  . 

BV0 

It  may  be  shown  that  u  >  0  everywhere,  as  is  required  of  equation  29.  First,  we  note  that 


(31) 


u(a  0)  —  Mo  — 


g sec  «o 

im 


(32) 


which  is  positive.  Then,  we  determine  that 

du 

da 


=  —2  sec3  a  , 


(33) 


which  is  everywhere  negative.  But  since  da  is  negative  along  the  trajectory,  it  implies  that  du  is 
positive,  and  that  the  initially  positive  value  u  =  uq  will  amplify  along  the  trajectory. 


7.  Ancillary  Velocity  Variables 


From  the  initial  integration,  resulting  in  equation  29  expressing  a  as  a  function  of  a,  a  number  of 
velocity  variables  may  be  immediately  established.  First,  from  equation  11,  the  projectile  speed 
may  be  expressed  in  terms  of  a  as 


s  = 


j~g~  sec  a 
V  B  y  /u 


(34) 


The  Cartesian  velocities  are  also  accessible  knowing,  along  the  trajectory,  that  dx  =  ds  cos  a  and 
dy  =  ds  sin  a.  It  immediately  follows  that 


x  = 


(35) 


and 


tan  a 


(36) 


For  completeness  sake,  we  present  the  time  rate  of  change  of  the  intermediate  variable  u,  given  as 


du 

u  =  a  —  = 
da 


2 \fBg-  \[u 


sec  a 


(37) 


At  this  stage,  the  specification  of  an  intermediate  angle  (a)  along  the  trajectory  of  the  projectile 
flight  is  sufficient  to  yield  the  associated  horizontal  and  vertical  velocities,  as  well  as  the  projectile 
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speed  along  the  trajectory.  Further,  the  time  rate  of  change  of  trajectory  angle  is  determined,  as 
well,  by  equation  29. 

The  one  unfortunate  result  is  that,  since  equation  28  cannot  be  algebraically  inverted  to  yield 
a(u),  one  is  not  able  to  determine  (in  closed  form)  the  associated  trajectory  angle  by  first 
specifying  a  velocity  along  the  trajectory. 


8.  Integrating  for  the  Trajectory  Pathlength 


According  to  the  chain  rule,  one  may  express  ds/du  as  the  ratio  of  s  and  u: 

ds  s  1 

du  u  2  Bu 


(38) 


Separating  the  variables  gives 


Straightforward  integration  gives 


ds 


s  = 


du 
“  2Bu 

In  (u/uq) 


2  B 


(39) 

(40) 


where  u0  has  been  previously  defined  in  equation  32.  Equation  40  provides  the  pathlength  s 
along  the  trajectory,  as  a  function  of  u  and,  by  inference,  as  a  function  of  the  trajectory  angle  a. 


The  analytical  determination  of  s  allows  other  useful  relations  to  be  algebraically  constructed,  for 
example, 


x 

T0 


=  e 


-Bs 


(41) 


9.  Known  and  Unknown 


We  have  demonstrated  a  solution  to  the  general  trajectory  problem  of  a  ballistic  projectile  subject 
to  aerodynamic  drag  in  a  gravitational  field.  The  problem  was  cast  in  curvilinear  coordinates  in 
order  to  obtain  the  analytical  solution.  As  a  result  of  the  two  integrations  of  the  second-order 
governing  equation,  we  were  able  to  obtain  the  trajectory  pathlength  s  as  a  function  of  the 
trajectory  angle  a. 
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Also  obtained  were  the  Cartesian  and  pathlength  velocity  components  (. x ,  y,  and  s,  respectively) 
as  well  as  the  time  rate  of  change  of  trajectory  angle  (a),  all  of  which  are  expressed  in  terms  of 
the  current  trajectory  angle  a. 


Unfortunately,  there  are  also  quantities,  which  would  be  very  useful  to  know,  for  which  the 
current  solution  does  not  provide.  The  intermediate  variable  u(a),  essential  to  the  problem 
solution,  may  not  be  analytically  inverted  for  a(u),  even  though  an  excellent  approximation  to  the 
inversion  exists*  in 


a(u)  ~  ±  arctan  ^ 


-1  +  \J\  +  0.3(C  —  u)‘ 
U6 


(42) 


which  is  within  0.37%  of  the  exact  solution  over  the  full  domain.  The  negative  result  is  used 
when  ( C  —  u)  is  negative.  However,  when  dealing  with  exact  results,  the  independent  variable  of 
the  problem  must  remain  a,  such  that  a  cannot  be  analytically  determined  as  a  function  of  any 
other  variable,  such  as  pathlength  s. 


Secondly,  while  velocity  components  were  obtained,  to  include  Cartesian  velocity  components  x 
and  y,  respectively,  these  are  not  analytically  integrable  using  the  current  approach.  Thus,  the 
solution  does  not  provide  the  Cartesian  coordinates  of  the  trajectory  (x  and  y),  as  an  analytical 
function  of  any  of  the  variables  in  the  problem.  However,  for  the  special  case  where  the 
trajectory  angle  remains  flat  throughout  ( a  ~  0),  the  horizontal  position  component  x  is 
adequately  described  by  the  trajectory  pathlength  s.  Also,  integrated  as  well  as  approximate 
analytical  solutions  in  Cartesian  coordinates  are  offered  later  in  this  report. 


Finally,  the  integration  of  a  over  time  could  not  be  analytically  achieved.  Thus,  while  many 
variables  are  known  in  terms  of  the  trajectory  parameters,  the  time  required  to  reach  various 
points  on  the  trajectory  is  not  analytically  known  from  the  current  approach,  even  though 
equation  29  may  be  separated,  with  the  time  being  precisely  expressed  in  terms  of  a  as 


t  = 


1  rQ  sec2  a 
\J  Bg  Ja0  \fu 


da 


(43) 


However,  for  the  special  case  where  the  trajectory  angle  remains  flat  (a  ~  0),  one  may  transform 
time  to  a  new  variable  r,  defined  by  the  transform  dr/dt  =  sec  a,  such  that  r  =  0  when  t  =  0. 
The  variable  r  may  be  explicitly  solved  and,  for  this  special  case  where  a  ~  0,  one  finds  i  «  r 
and  gets 


t 


<A~0 


T 


VU  -  y/U^ 

VBg 


(44) 


*Equation  42  was  developed  by  determining  that  (C  —  u)  =  tan  a  sec  a  +  ln(tan  a  +  sec  a)  was  extremely  well 
approximated  by  (C  —  u)  «  2  tan  a(Vl  +  0.3  tan2  a).  This  approximation  was  inverted  to  yield  equation  42. 
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10.  Expedited  Numerical  Integration  for  Cartesian  Trajectory 


Because  the  trajectory  s  is  known  in  terms  of  a,  an  expedited  numerical  integration  may  be  set  up 
that  provides  the  time  history  of  the  trajectory  in  Cartesian  coordinates.  Instead  of  timesteps,  the 
method  takes  small  negative  ct: -steps  from  the  original  trajectory  a0.  At  launch,  the  values  of  a, 
u,  and  s  are  a0,  u0,  and  V0,  respectively.  After  a  is  decremented  during  an  o-stcp,  the  updated 
values  for  u  and  s  are  analytically  obtained  from  equations  28  and  34,  respectively.  The  exact 
trajectory  pathlength  increment  associated  with  the  o-stcp  is  analytically  obtained  as 


ds  =  s  —  s 


In  (u/u  ) 
2 B 


(45) 


where  the  “minus”  superscript  denotes  values  at  the  prior  o-stcp.  Directional  components  of  this 
increment  are  used  to  update  the  Cartesian  (x,  y )  coordinates  from  their  initial  values  of  (0,  II). 


x  =  x  +  ds  ■  cos  a 


(46) 


y  =  y  +  ds  ■  sin  a  ,  (47) 

where  a  =  (a  +  a~)/2.  The  time  associated  with  a  is  obtained  from  ds  and  s  as  follows: 

t  =  t~  +  ds/s  ,  (48) 

where  s  =  (s  +  s~)/2.  In  this  simple  fashion,  with  successive  o-stcps,  the  time  history  of  the 
trajectory  may  be  obtained  in  Cartesian  coordinates.  A  FORTRAN  listing  that  implements  this 
algorithm  is  given  in  table  1. 

The  results  of  this  algorithm  are  demonstrated  in  figures  2-4.  The  simulation  was  for  the 
trajectory  of  a  50  mm-long  copper  (pp  =  8900  kg/m3)  projectile  with  a  drag  coefficient  of  1.0, 
launched  from  an  initial  altitude  of  H  =  1  m  at  a  velocity  of  2000  m/s  and  a  launch  angle  of  25°. 
For  these  initial  conditions,  the  value  for  the  parameter  B  is  0.00145  m_1.  In  these  figures,  a 
symbol  is  plotted  for  every  30  m  of  pathlength  traveled.  Recall  that  only  x,  y,  and  t  are 
numerically  integrated  quantities.  All  other  plotted  values  (including  s,  s,  a,  and  da  Ids)  are 
analytically  defined  in  closed  form.  Of  course,  the  same  result  could  have  been  achieved  through 
a  conventional  discretization  and  integration  of  the  governing  equations,  equations  5  and  6. 
However,  the  advantage  of  the  algorithm  given  in  table  1  is  computational  speed  and  algorithmic 
simplicity. 
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Table  1.  FORTRAN  code  listing 


C  TRAJECTORY  MODEL,  WITH  NUMERICAL  INTEGRATION  FOR  CARTESIAN  QUANTITIES 

c 

implicit  none 
integer  n 

double  precision  B,  VO,  g,  a,  u,  aO,  H,  da,  C,  uO,  PI,  af,  ds, 

&  x,  y,  t,  uget,  uf,  sdot,  s,  dsprint,  sdotf,  adot 

data  g,  da,  n,  dsprint  /9.806,  -0.001,  1,  30./ 

PARAMETER  (PI  =  3.1415926) 

print  *,  'Enter  B,  Vo,  alphaO,  H:  ' 
read  (*,*)  B,  VO,  aO,  H 

write (7, '  ( ' ' B,  VO,  aO,  H:  ' ' , lp4ell . 4 , /) ' )  B,  VO,  aO,  H 
write  (7, ' '  V  alpha  s  x' ' , 

&  '  '  y  t  da/dt' ' ) ' ) 

aO  =  aO  *  PI/180, 
da  =  da  *  PI/180. 

uO  =  g  /  (B  *  V0**2  *  (cos (aO) ) **2) 

C  =  -uget(a0,  -uO) 

PRINT  *,  'uO,  C:  ' ,  uO, C 

C  INITIALIZE 

a  =  aO 
x  =  0. 

y  =  H 

s  =  0 . 
sdot  =  VO 
t  =  0. 

u  =  uget (a,  C) 

adot  =  -sqrt (B*g*u) * (cos (a) ) **2 

write  (7,' (7fl0.3)')  sdot,  a  *180. /PI,  s,  x,  y,  t,  adot 

C  ALP HA- STEP  LOOP : 

do  while  (y  .ge.  0.) 

C  SAVE  PRIOR  ALPHA-STEP  QUANTITIES;  THEN  STEP  ALPHA 

uf  =  u 
sdotf  =  sdot 
a  =  a  +  da 

C  GET  CURRENT  ALPHA-STEP  QUANTITIES  ANALYTICALLY 

u  =  uget (a,  C) 
s  =  log(u/u0)  /  (2.*B) 
sdot  =  sqrt (g/ (B*u) )  /  cos (a) 
ds  =  log(u/uf)  /  (2 . *B) 
adot  =  -sqrt (B*g*u) * (cos (a) ) **2 

C  HERE  ARE  THE  INTEGRATED  QUANTITIES 

x  =  x  +  ds  *  cos  (a  -  da/2.) 
y  =  y  +  ds  *  sin (a  -  da/2.) 
t  =  t  +  2.*ds/(sdot  +  sdotf) 

C  CHECK  IF  TIME  TO  PRINT 

if  (s  .ge.  float (n) *dsprint)  then 

write  (7,'(7fl0.3)')  sdot,  a  *  180. /PI,  s,  x,  y,  t,  adot 
n  =  n  +  1 
end  if 

if  (abs (a)  . le .  -da/2.) 

&  write  (7,'(7fl0.3)')  sdot,  a  *  180. /PI,  s,  x,  y,  t,  adot 
end  do 

write  (7,'(7fl0.3)')  sdot,  a  *  180. /PI,  s,  x,  y,  t,  adot 

stop 

end 

c* ************************************* 

double  precision  function  uget (a,  C) 

implicit  none 

double  precision  a,  C 

uget  =  C  —  tan (a) /cos (a)  -  log (tan (a)  +  l./cos(a)) 

return 

end 


11 


Figure  2.  Trajectory  (y)  and  velocity  ( ds/dt )  for  50  mm-long  copper  (pp  =  8900 
kg/m3)  projectile  with  Cd  =  1-0  (B  =  0.00145/m)  launched  at 
2000  m/s  at  an  angle  of  ao  =  25°  from  an  initial  altitude  of  H  =  1  m. 
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11.  Special  Case  Solutions  and  Approximations 


11.1  Approximate  Cartesian  Solution  for  Moderate  Angle-of-attack  Trajectories 

While  the  exact  trajectory  solution  has  been  presented  in  curvilinear  coordinates  (s,  a),  it  was 
made  clear  in  section  9  that  a  corresponding  analytical  solution  was  not  achieved  in  Cartesian 
coordinates  ( x ,  y ).  Nonetheless,  an  excellent  approximation  can  be  derived,  which  is  valid  for 
trajectories  through  moderate  angles  of  attack. 


(49) 

(50) 

(51) 


While  we  have  been  unable  to  solve  these  integrals  when  u  is  accurately  defined  by  equation  28,  it 
is  possible  to  develop  an  approximate  solution  by  arbitrarily  redefining  an  approximation  to  u 
(call  it  u)  that  is  both  integrable  in  the  context  of  equations  49-5 1  as  well  as  approximately  equal 
to  equation  28  over  a  substantial  range  of  the  functional  domain. 


While  there  are  many  non-logarithmic  functions  that  can  be  employed  to  very  accurately 
approximate  equation  28,  their  substitution  into  equations  49-51  must  result  in  analytically 
integrable  functions.  Thus,  the  choices  are  quite  limited,  and  the  function  we  propose  is 


u  =  C  —  2  tan  a 


There  are  two  facets  to  this  approximation  to  equation  28.  One  is  the  use  of  the  constant  C  as 
opposed  to  C.  The  other  is  the  use  of  2  tan  a  in  place  of  tan  a  sec  a  +  ln(tan  a  +  sec  a).  As  to 
C,  its  specification,  in  order  to  rigorously  satisfy  the  initial  t  =  0  condition,  should  be 


-  g  sec  a0 

C  (nominal)  =  uq  +  2  tan  ao  =  — - b  2  tan  a0 

BV0 
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However,  since  this  development  is  part  of  an  approximate  solution,  we  are  willing  to  relax  rigor 
to  achieve  a  better  fit,  as  shown  later. 

For  the  second  facet  of  the  approximation,  we  see  in  figure  5  that  the  approximation  given  by 
2  tan  a  is  quite  good  up  to  ±20°,  and  maybe  acceptably  good  to  a  range  of  ±40°.  This  range  of 
applicability  would  cover  quite  a  span  of  useful  trajectories. 
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Figure  5.  A  comparison  of  tan  a  sec  a  +  ln(tan  a  +  sec  a ) 
to  the  approximation  given  by  2  tan  a. 


Turning  to  equations  49-51,  and  employing  the  approximation  to  u  given  by  equation  52,  one 
may  use  the  substitution  w  =  tan  a  to  reexpress  these  governing  equations  as 


(54) 


and 


dx 


1  dw 
B  C-  2w 


dy 


1  w  dw 
B  C-2w 


(55) 

(56) 
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These  equations  may  be  integrated  directly  to  give  the  following: 


t  = 


VB~g 


C  —  2w  —  \  C  —  2 w, 


(57) 


X  ~  2B  hl 


and 


y 


H  = 


w 


'  C  —  2w 
v  C  —  2  w0  , 


-W0  C 
2B~  +  lBln 


C  —  2w 


KC  —  2wq  J 

where  w  and  w0  are,  respectively,  tan  a  and  tan  a0;  and  II  is  the  altitude  of  the  launch  point. 
The  equations  for  x  and  y  may  be  alternately  expressed,  eliminating  w,  as 


(58) 


(59) 


VBgt  \ 
C  —  2wq  ) 


(60) 


y  = 


Cx 

2 


gt 2 

4 


+  H 


(61) 


We  note  that  these  alternate  representations  mimic  the  functional  dependencies  derived  in  an 
earlier  shallow-trajectory  model  (5),  though  with  different  constants. 


These  approximate  solutions  were  compared  against  the  integrated  trajectory  of  equations  46-47. 
The  same  trajectory  problem  was  solved  as  presented  in  section  10,  except  for  varying  the  initial 
trajectory  angle,  a0.  On  the  basis  of  those  results,  the  value  of  C  constant  was  refined  from  that 
given  in  equation  53  to  a  value  of 


C  =  «o(l  +  2.5  tan2  ck0)  +  2tana0  , 


(62) 


in  an  effort  to  extend  the  applicability  of  the  fit  to  larger  angles  of  attack.  The  effect  of  modifying 
C  from  its  nominal  value  is  to  introduce  a  less  accurate  value  of  a  at  t  =  0,  with  the  goal  of 
improving  the  late-time  trajectory  approximation.^ 

With  the  approximations  developed  in  equations  57-59  and  the  fit  to  C  given  by  equation  62,  the 
trajectories  are  compared  in  figures  6-7.  We  see  that,  indeed,  the  approximation  is  excellent  for 
initial  trajectories  at  or  below  20°,  with  trajectory  range  errors  well  below  1%,  and  altitude  errors 
below  1.5%.  The  flight  duration  approximation  was  in  error  by  under  4%,  as  well.  As  the 

^Equation  62  was  refined  for  this  particular  test  case.  A  more  general  form,  which  tries  to  account  for  vari¬ 
abilities  in  B  and  V0,  in  addition  to  a0,  is  given  by  C  =  Uq { 1  +  9/8(T4o/V6o)3/ (1  +  (fo/Voo)3)  [5.906  + 
1.21og10(B/B0)]  tan2  ct0}  +  2tana0,  where  B0  =  1  m-1  and  Too  =  1000  m/s. 
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Figure  6.  A  comparison  of  numerically  integrated  (black) 
and  approximated  (red)  trajectories  for  initial 
trajectory  angles  below  20°. 


Figure  7.  A  comparison  of  numerically  integrated  (black) 
and  approximated  (red)  trajectories  for  initial 
trajectory  angles  between  25°  and  40°. 
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launch  angle  is  increased  beyond  20°,  the  approximation  is  less  accurate  at  predicting  total  range, 
even  as  the  approximation  of  maximum  altitude  remains  of  good  quality. 

But  even  so,  at  a  30°  launch  angle,  the  integrated  maximum  altitude  is  941  m,  whereas  the 
approximated  maximum  altitude  is  928  m,  a  1.4%  difference.  The  maximum  horizontal  range  of 
that  condition  is  integrated  as  2578  m,  and  approximated  as  2646  m,  a  mere  2.6%  difference. 

The  duration  of  the  flight  is  26.1  s  when  integrated,  compared  with  24.5  s  when  approximated,  a 
6.2%  difference. 


One  notes  that  the  approximation  loses  its  greatest  accuracy  during  the  terminal  phase  of  the 
trajectory.  This  can  be  explained  through  an  examination  of  the  original  approximation  found  in 
figure  5.  While  a0  starts  out  in  a  reasonable  part  of  the  fit,  during  the  trajectory,  a  is 
monotonically  decreasing.  After  the  maximal  altitude  is  reached  at  a  =  0,  the  value  of  a 
proceeds  rapidly  into  the  negative  domain  and,  in  the  terminal  phase  of  the  trajectory,  is  reaching 
values  in  the  vicinity  of  —80°,  where  the  approximation  offered  by  2  tan  a  is  particularly  poor. 
However,  even  with  these  drawbacks  at  higher  angles  of  attack,  the  results  of  the  current 
Cartesian  approximation  (given  by  equations  57-59)  were  confirmed  as  more  accurate  than  the 
shallow-angle  trajectory  approximation  given  in  reference  5. 

11.2  Maximum  Pathlength  from  Launch  to  Apogee 


While  it  is  true  that  the  exact  trajectory  solution  was  not  expressed  in  Cartesian  coordinates,  one 
may  nonetheless  employ  the  curvilinear  solution  to  glean  interesting  facets  of  the  trajectory.  For 
example,  the  apogee  is  the  point  in  the  trajectory  of  maximum  altitude  (where  a  =  0).  Since  the 
value  of  u  equals  C  at  the  apogee  (equation  28),  it  immediately  follows  from  equation  35  that  the 
velocity  at  the  apogee  is  ^Jg/BC. 

Whereas  such  details  may  be  obtained  through  direct  algebraic  substitution,  more  complex 
information  may  also  be  derived  from  the  solution.  For  example,  the  analytical  solution  may  be 
employed  to  determine  the  launch  angle  needed  to  maximum  the  pathlength  from  launch  point  to 
apogee,  as  we  demonstrate  here. 


The  pathlength  traversed,  from  launch  to  apogee,  call  it  sap,  is  obtained  from  substituting 
equations  3 1  and  32  into  equation  40: 


In  (C/u0)  1 

s„„  =  - =  —  in 

p  2  B  2  B 


1  + 


BVn  /tan «0 


H - ; 

g  \  sec  ao  sec^  ao 


ln(tan  a0  +  sec  a0) 


(63) 


To  find  the  launch  angle  a0  that  maximizes  sap,  one  merely  need  take  dsapj daM  and  set  the  result 
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to  zero.  That  process  leads  to  the  result  that  the  maximum  apogee  pathlength  is  obtained  for  that 
value  of  oto,  which  satisfies  the  following  relationship: 

ln(tanao  +  secao)  =  csccto  [  for  max(sap)  ]  .  (64) 


This  relationship  is  satisfied  when  a0  =  56.4656  . . ,°,  and  is  a  very  surprising  relationship.  It 
indicates  that,  regardless  of  the  drag  coefficient  embodied  in  B,  the  initial  launch  velocity  Vo,  and 
the  gravitational  constant  g,  the  launch  angle  that  maximizes  the  apogee  pathlength  is  always 
56.5°!  Indeed,  through  equation  63,  the  pathlength  itself  depends  upon  B,  V0,  and  g,  but  the 
launch  angle  to  maximize  that  apogee  pathlength  does  not.  This  result  has  been  computationally 
verified  for  trajectories  under  the  influence  of  aerodynamic  drag  and  analytically  verified  for 
drag-free  (parabolic)  trajectories. 

11.3  Terminal  Trajectory 

When  a  trajectory  is  in  its  terminal  stages,  given  enough  altitude,  it  will  approach  the  terminal 
velocity,  which  is  given  by  \J g/B.  An  examination  of  equation  34  for  the  projectile  velocity,  s, 
reveals  that,  in  the  terminal  phase,  y/u  — >  sec  a,  in  order  for  the  velocity  to  approach  its  terminal 
value. 


Therefore,  if  we  substitute  sec  a  for  the  appearance  of  y/u,  thereby  forcing  the  projectile  velocity 
to  remain  fixed  at  the  terminal  velocity,  we  should  be  able  to  leam  something  of  the  terminal 
trajectory.  Let  us  use  the  t  subscript  to  define  variables  at  the  onset  of  the  terminal  phase. 


First,  consider  equation  29,  which  defines  a.  Once  yfu  is  replaced  with  sec  a,  the  equation  may 
be  separated  to  solve  for  dt  and  then  integrated: 


1 

y/Bg 


sec  a  da 


resulting  in 


t  —  tf  = 


In 


sec  at  +  tan  at 
sec  a  +  tan  a 


Recall  that,  in  the  terminal  phase,  a  will  be  negative,  and  tending  towards  —90°. 


(65) 


(66) 


Turning  to  the  x  coordinate  of  the  trajectory,  the  same  approximation  may  be  made  to  equation  35 
to  obtain 


(67) 
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leading  to 


Similarly  for  y. 


leading  to 


x  -  xt  =  -  —  (a  -  at) 

ID 


ry  l  ra 

/  dy  =  ——  tan  a  da  , 
lyt  -D  J at 


1  /sec \ 

y  -  2/t  =  n  ln  - 

B  \  sec  a  j 


(68) 

(69) 

(70) 


These  results  bear  some  similarity  to  the  1920  report  by  Wilson  (1)  on  bomb  trajectories,  in  which 


the  y  trajectory  is  approximated  in  various  ways,  the  sixth  such  approximation  being  expressed 
largely  in  terms  of  a  —  ln(sec  x)  term. 


To  see  how  this  function  of  terminal  trajectory  compares  to  the  exact  solution,  we  provide 
figure  8,  which  zooms  in  on  the  terminal  phase  of  the  trajectory  described  in  figure  2. 
Superimposed  on  the  graph  are  three  instances  of  the  approximated  terminal  trajectory  implicitly 
given  by  equations  68  and  70.  In  these  three  terminal  trajectories,  the  starting  point  associated 
with  (xt,  ry)  corresponds,  respectively,  to  a,  values  of  0°,  —30°,  and  —60°.  Obviously,  the  closer 
at  is  taken  to  —90°,  the  more  accurate  the  estimate  of  the  terminal  trajectory  will  be. 
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Figure  8.  Estimated  terminal  trajectories  from  three 
stalling  points  on  the  actual  trajectory, 
corresponding  to  at  equal  to  0°,  —30°,  and 
—60°,  respectively. 
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12.  Conclusions 


The  solution  to  the  problem  of  low-altitude  ballistic  trajectory,  under  the  influence  of  a  uniform 
gravitational  field  and  aerodynamic  drag,  has  been  presented  here.  By  “ballistic,”  we  are  limiting 
the  analysis  to  a  projectile  that  possesses  neither  propulsion  nor  lift.  The  solution  was  achieved 
through  a  transformation  of  the  governing  equations  into  a  curvilinear  coordinate  system.  The 
solution  provides  for  the  pathlength  of  the  trajectory  in  terms  of  the  trajectory  angle.  Also 
available  from  the  solution  are  the  velocities  along  the  pathlength;  the  vertical  and  horizontal 
velocity  components;  and  the  time  rate  of  change  of  the  trajectory  angle.  In  all  cases,  these 
variables  are  expressed  in  terms  of  the  current  trajectory  angle. 

Unfortunately,  there  are  useful  functions  of  interest  that  are  not  analytically  accessible  from  the 
solution.  These  include  the  Cartesian  (x,  y )  coordinates  of  the  trajectory  and  the  time  required  to 
reach  various  stages  of  the  trajectory.  Also,  the  solution  cannot  be  mathematically  inverted.  As 
a  result,  the  trajectory  angle  must  remain  the  independent  variable  of  the  problem,  and  cannot  be 
analytically  expressed  in  terms  of  other  quantities  along  the  trajectory. 

To  increase  the  utility  of  the  solution,  when  a  Cartesian  solution  is  essential,  an  integration 
algorithm  was  presented,  which  makes  use  of  the  exact  solution,  in  order  to  expedite  the 
calculation  of  the  Cartesian  time  history  of  the  trajectory.  In  the  algorithm  presented,  the 
discretized  integration  variable  is  not  time,  as  would  traditionally  be  the  case  when  employing  a 
brute-force  approach  to  the  problem,  but  rather  the  trajectory  angle.  The  algorithm  was 
demonstrated  on  a  sample  problem,  and  the  results  were  graphically  presented. 

Several  special-case  solutions  were  offered,  including  an  analytical  approximation  to  the 
Cartesian  solution,  quite  accurate  up  to  moderate  angle-of-attack  trajectories  (less  than  20°, 
approximately),  and  reasonably  accurate  beyond  that  (e.g.,  up  to  30°).  Further,  an  unexpected 
analytical  result  was  developed,  which  shows  that  the  trajectory  that  maximizes  the 
apogee-pathlength  ( i.e .,  pathlength  from  launch  to  apogee)  is  56.5°,  regardless  of  launch  velocity, 
drag  coefficient  (including  drag-free  parabolic  trajectories),  and  gravitational  constant. 
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